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o\ : 

^\ i A single spin flip stochastic energy conserving dynamics for the XY model is considered. We 

study the ergodicity properties of the dynamics. It is shown that phase space trajectories densely 
» I ' fill the geometrically connected parts of the energy surface. We also show that while the dynamics is 

f> |_ discrete and the phase point jumps around, it cannot make transitions between closed disconnected 

parts of the energy surface. Thus the number of distinct sectors depends on the number of geo- 
metrically disconnected parts of the energy surface. Information on the connectivity of the surfaces 
is obtained by studying the critical points of the energy function. We study in detail the case of 
two spins and find that the number of sectors can be either one or two, depending on the external 
fields and the energy. For a periodic lattice in d dimensions, we find regions in phase space where 
the dynamics is non-ergodic and obtain a lower bound on the number of disconnected sectors. We 
O ' provide some numerical evidence which suggests that such regions might be of small measure so 

that the dynamics is effectively ergodic. 
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I. INTRODUCTION 



A much studied problem in non-equilibrium statistical mechanics is the problem of heat conduction. Here one 
is interested in obtaining the steady state of a system whose two ends are kept at different temperatures. Some 
of the important problems are to see if local thermal equilibrium is achieved in the steady state, to calculate the 
£j ' temperature profile inside the system and to check the validity of linear response laws relating heat current and 
temperature gradient (Fourier Heat Law). In this context models with both deterministic and stochastic dynamics 
have been studied. Rieder et al exactly solved the problem of heat conduction in harmonic chains M. They found a 
, constant temperature distribution inside the system and a heat current that was independent of system size. Since 
then, various authors have studied the effect of impurity and nonlincarity in one and two dimensional systems of 
interacting atoms evolving via deterministic Newtonian dynamics . Numerical studies show that with some kinds 
of interaction one gets realistic temperature profiles and heat currents proportional to inverse of system size Q]. For 
the Lorenz gas in the Boltzmann-Grad limit, it has been proved that the linear Boltzmann equation holds and in 
this case too one gets the Fourier Law ||. There are however not too many studies on heat conduction in stochastic 
models. Connor and Lebowitz studied the random reflection model in which particles moving on a line are reflected 
with some probability from randomly spaced scatterers j|. Streater has studied heat conduction in a lattice model 
with finite degrees of freedom ||. 

In this paper we are interested in stochastic lattice models of heat conduction. Specifically we look at classical spin 
systems with some kind of imposed dynamics. In studying spin systems one typically uses dynamics such as Glauber 
where each point of the lattice is kept in contact with heat baths at given temperature. In this case the spin flip 
rates are assumed to satisfy detailed balance condition, with the rates depending on the heat bath temperature. For 
equilibrium studies the same temperature is imposed at all points and the steady state is the canonical distribution. 
However, in studying heat conduction, the temperature distribution in the steady state is not known beforehand and 
should come out of the dynamics. Since the temperature is space dependent and not known in the beginning, one 
cannot specify the relaxation rates for energy raising and energy lowering transitions. The only simple choice is to 



5— I 

make such rates zero. Thus we use dynamics where there is no exchange of heat with external reservoirs (except at the 
two ends) and there is local conservation of energy. However, in this case the dynamics becomes highly constrained, 
and it is not clear that it wont get stuck. A requirement of correct phenomenological description is that the system 
should relax to a state of local thermal equilibrium. In order that this occurs, we expect that the isolated system 
should reach thermal equilibrium, that is the dynamics should be ergodic and lead to a microcanonical distribution. 
In this paper we address this problem only. We study heat conduction in our model in a future paper. 

It is easy to show that the Ising model with single spin flip energy conserving dynamics is non-ergodic in dimensions 
higher than one. For example, suppose we take a (i-dimensional hypercubic lattice of N sites, and consider the spin 
configuration of one single up-spin in a sea of down-spins. If we take periodic boundary conditions then there are 
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N degenerate energy configurations corresponding to the N places where the single up-spin can be located. The 
dynamics does not connect these states. For other energies too one can show that there is no ergodicity. In order to 
do microcanonical simulations for the Ising model it is thus necessary to couple it to other degrees of freedom with 
which it can exchange energy. This has been studied by Bhanot et al [Q. Heat conduction in the Ising model (with 
additional degrees of freedom) was studied by Creutz ||. 

One can ask what happens if we replace the Ising spins by continuous spins. Microcanonical simulations of the 
XY model have earlier been done by introducing additional degrees of freedom. This include the molecular dynamics 
study by Kogut and Polonyi who introduce momentum variables and use a continuous time Hamiltonian dyanmics 
H . A microcanonical simulation using the Creutz algorithm was done by Ota et al . 

In this paper we consider an energy conserving single spin flip dynamics for the XY model, with no additional 
degrees of freedom, and look at its ergodicity properties. We first study the instructive case of two spins interacting 
with each other and with external fields. This case is still non-trivial but can be worked out in full detail. We show 
that for generic values of external fields, the phase space trajectories (note that the trajectory is not continuous and is 
formed by the random iterates of a phase point) densely fill the geometrically connected parts of the constant energy 
contours. Thus whenever the energy contours are fully connected wc have ergodicity. On the other hand, even though 
the dynamics is discontinuous, we show that jumps between closed disconnected parts of the energy contours are not 
allowed. Thus in such cases we have non-ergodic behavior and we have as many sectors as disconnected energy parts, 
within each of which the dynamics is ergodic. Thus ergodicity is related to nature of the constant energy contours. 
Information regarding the nature of the energy contours can be obtained from a study of the critical points. These 
are points at which the derivative of the energy w.r.t all angles vanish. For instance, if there are two or more stable 
critical points (these correspond to fixed points of the dynamics) with the same energy values then we will have closed 
energy contours around each of them, which are disconnected. We find that the nature of energy contours depend on 
the external fields. For some fields there are only single sectors at all energies while for other field values there are 
one or two sectors depending on the energy. 

Next we study the case of N > 2 spins. Using the two spin results, we prove that even for the general case of N 
interacting spins, the phase space trajectory densely covers the geometrically connected parts of the constant energy 
surface. Thus ergodicity is related to the connectivity properties of the energy surfaces. As in the two-spin case, some 
information can be obtained by looking at the critical points. As an illustrative example, we study a system of three 
spins interacting with each other and with external fields. We show that the energy surfaces can be constructed by 
studying the critical points. In this case all energy surfaces are connected and phase space trajectories densely fill the 
surfaces. Finally we study spins on hypercubic lattices in d dimensions. In this case we are not able to obtain a full 
characterization of the energy surfaces. However we explicitly construct a large number of spin configurations which 
correspond to critical points, both stable and hyperbolic, and which span the entire energy spectrum of the system. 
The presence of degenerate stable points implies disconnected energy surfaces around the neighborhood of these 
points. Thus in such cases we get non-ergodic behavior. However we find some numerical evidence which suggests 
that such regions where the dynamics is non-ergodic are of small measure. For example, we find from simulations 
using the energy conserving dynamics that long time averages of correlation functions are the same as that obtained 
from a canonical simulation (which can be shown to be ergodic) whenever the initial conditions are chosen randomly. 
For initial conditions close to stable points the correlation functions can be calculated exactly and differ from the 
canonical averages. On the other hand if we start with initial conditions close to hyperbolic points we get from 
simulations results close to the canonical ones. 

The paper is organized as follows. In section II wc define the model and the dynamics. In section III the two spin 
case is first studied. The general case of N > 2 spin case is studied in section IV. 



In the XY model each spin is a two-component unit vector specified by the angle 9 that it makes with some fixed 
axis. Equivalently every spin can be denoted by the complex number z — e 10 . The interaction energy for spins on a 
lattice is given by 



II. DEFINITION OF MODEL 




(1) 
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where the sum is over all nearest neighbors. The dynamics is the following: during an infinitesimal time dt each 
spin undergoes an energy conserving flip with a probability dt. For discrete time evolution, as would be required for 
a computer implementation of the dynamics, this is equivalent to choosing a spin randomly at every time step and 
flipping it. Let h = J2ien n z « ^ e tne net magnetic field at a site due to all the other nearest neighbor sites. Then the 
spin flip is given by the map 

z' = hz*/h*. (2) 

which corresponds to a reflection of the spin z about the direction h. 

For ergodic behavior it is necessary that starting from a generic initial condition it should be possible to access all 
other points on the constant energy surface. The detailed balance condition then ensures that a uniform probability 
distribution over the energy surface is the stationary distribution. 

We note that a generalization of the above dynamics to the n-vector model can be made. For the Ising model 
(n = 1) a flip can only occur if a spin has equal number of up and down neighbors. For the XY model (n = 2) there 
is always a unique flipped state. For higher n the spin can flip to a continuum of values. 



III. ERGODICITY IN TWO SPIN SYSTEMS 



Consider two spins z\ and z 2 coupled to each other by the xy-interaction and also interacting with constant external 
fields of magnitude h\ and h 2 and directions given by the angles fa and fa. The fields h\ and h 2 could be the effect 
of other spins held fixed in a bigger system. The energy of the system is given by 

Ti = —hicos(9i — cf)i) — h2Cos(9 2 — fa) — cos(8i — 9 2 ). 

Here 9\ and 9 2 are the angles of the two spins. We can choose the reference axis such that fa — and fa = <j>. If 
we are interested only in finding the region that is accessed by the system , then for the two spin case, we may as 
well flip the spins successively. This is sufficient since flipping a spin twice successively leaves it unchanged. Thus we 
study the deterministic maps 



hi + z 2 , _ bgj* 

h 1 + l/zi /Zl > Z *- h2e -i<t> + 



We first note the following points. In general any constant 9i line will intersect the constant energy contour given 
by H = e at two points only. These correspond to the two equal energy configurations of a spin in a net magnetic 
field. The energy contours can thus be of the forms shown in fig. 1. In (a) the contours are connected while in (b) 
and (c) they are disconnected. However it is clear that the dynamics allows transitions between the disconnected 
parts in (b) but not so in case (c). This is because, at every step of the dynamics, only one spin can change. But 
moving between the disconnected parts of the energy contours in (c) requires two spins to change simultaneously. 
Thus whenever there are closed disconnected contours, that are non-spanning, there can be no transitions from one 
to the other. 

We now show that for generic fields and initial conditions the trajectories are non-periodic and densely fill the 
constant energy contours. The proof consists of the following observations: 

(i) Considered as a map in the 0\, 9 2 plane, (3) is area preserving, i.e it has jacobian of modulus one. This is easily 
proved by writing eqns. (3) in the form 9[ = —6\ + f(9 2 ),9' 2 = —9 2 + g(9[), where / and g are given functions. Thus 
the map can have no attractive periodic points. 

(ii) The map (3) is equivalent to a smooth invertible map on the circle. Since there can be no attracting set of 
points, two kinds of behavior of the trajectories is possible. Either all trajectories are periodic with the same period 
or else trajectories are aperiodic and densely fill the circle. 

(iii) For periodic orbits of order n we must have z[ n ^ = z\ and z^ — z 2 . For generic values of the fields, these 
equations will have a finite number of solutions, and therefore it follows that all points cannot belong to periodic orbits. 
Hence, because of (ii), all trajectories are aperiodic and fill the energy contours densely. It is however possible that 
at some special value of the fields the above eqns for periodic orbits are satisfied identically for some n. In that case 
all trajectories are periodic orbits of order n. Thus, as an example, for h\ = h 2 = 1, the equations z^ = zi.z^ = z 2 
are identically satisfied. All points are then periodic of order 3 and there are is no ergodicity. However if we change 
h 2 by an infinitesimal amount then the trajectories become dense. This can be seen in fig. 2 where we plot phase 
space trajectories for different initial conditions for h\ = 1.0, h 2 = 1.001. 
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For the two spin case the nature of the constant energy contours can be completely determined by studying the 
critical points of the energy function H = —hi cos(#i) — h 2 cos(9 2 — <P) — cos($i — 9 2 ). The critical points are obtained 
by setting the first derivatives of H w.r.t 9i and 9 2 to zero. Thus we get: 

— = hi sin(0i) + sin(6»i - 9 2 ) = 

^=h 2 sin(9 2 - 4>) - sm{6i - 62) = 
oa 2 

(4) 

These equations can be written in the more transparent form: 

2 hi + z 2 

Zi = 



hi + l/z 2 
2 _ h 2 e^ + gi 

Z2 ~ h 2 e-^ + l/zr (5j 

This form also follows from eqns (3) and imply that at the critical points either the field at any site vanishes or 
is parallel to the spin at that site. If the field does not vanish at both the sites then the critical point will be a 
fixed point of the dynamics. The nature of the critical points can be found by finding the eigenvalues of the matrix 
Mki = d 2 TL/ dQkdQi evaluated at the critical points. We are interested in finding the phase diagram in the (hi,h 2 , 0)) 
space showing regions where the energy contours are connected and the regions where they are disconnected. As 
we shall show now this can be obtained by studying solutions of eqn. (5). First we note certain symmetries of eqn. 
(5). It is left unchanged by the following transformations: (l)ft-i — > —hi, z 2 — > — z 2 ; (2) h 2 — > — h 2 , Zi — > —z\; 
(3)zi -» 1/zi, z 2 -> l/z 2 , 4> -> -0; (4) -> 7r - 0, zi -> -1/zi, z 2 -> I/Z2; (5) hi -> /12, /12 -> zi e^/z 2 , 
z 2 — > e^/zi . It can also be shown that the equations determining the nature of the critical points are also invariant 
under the above transformations. Because of these symmetries it is sufficient to consider the subspace given by 
hi > h 2 > 0,0 < (f> < it/2. 

From eqn (5) we get the following equation for z\ 

h^e-^zf + (h\h\ + hl- hle'^zl + (h^e^ - 2/if/i 2 e"^)^ + {-1h\h\ - 2h\ + 2h\)z\ + 

(hjh^-^ - 2h\h 2 e^)zt + (h\h\ + h\ - hje^zi + hjh^ = (6) 

For every solution z\ of (6), Z2 is uniquely determined through the equation: 

h 2 z 1 (e i4 ' - e'^zl) 



z 2 



hi{z\ - \){h 2 e-^zi + 1) 



Thus in general we obtain six sets of solutions to (5). The solutions for which z\ and z 2 have modulus one correspond 
to the physically relevant solutions. 

Let us first consider the case = 0, for which (6) can be solved for z\. We obtain the following six solutions: 
(zi,z 2 ) = (l,l),(l,-l),(-l,l),(-l,-l),(e Wl ,e- W2 ) and (e"^ 1 ,^ 2 ) where cos(6»i) = {h\ - h\- h\h\)l2h\h 2 and 
cos(0 2 ) = {h\ — h\ — h\h\^)j2h\hi. We see that two of the critical points merge and disappear as the lines h 2 = 
hi/{h\ ± l),/i2 = hi/(l — hi) are crossed. In fig. 3 we plot the phase diagram which shows the nature of the critical 
points in different regions in the hi, h 2 plane. The important point to note is that in the region enclosed by the three 
curves, where we have three stable fixed points, we get disconnected energy contours. This is because of the following 
reason. Consider what happens as we move along some constant hi line as shown in fig. 3. At point A just below the 
phase boundary we have two stable critical points; (1, 1), which is a energy minima and (—1,1) which is a maxima. 
The other two points (—1,-1) and (1,-1) are hyperbolic points. As we cross the line to point B, (—1,1) becomes 
unstable and at the same time, two new stable points, degenerate in energy appear from this point. The separatrix 
passing through (—1, 1) encloses both the stable points. This at once implies that there will be disconnected closed 
energy contours of identical energies around the two stable points. In fact for all energy values greater than the 
separatrix energy, we get disconnected sectors, while for lower energies there are only single sectors. Similarly as we 
cross from the point D to E in fig. 3, the two degenerate stable maxima merge with the unstable point (—1, —1) to 
give rise to a single stable point. Hence at E and all points above this we have only connected energy contours for 
all energy values. To illustrate the above point we plot in fig. 4 the energy contours for field values corresponding 
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to the points A, B, C, D and E. The energy contours were obtained by starting from different initial conditions and 
evolving through the map (3). 

For 4> ^ 0, we can no longer solve eqn (6) explicitly. However we can still obtain the phase diagram. In this case 
for fixed h\ suppose we increase the h 2 value starting from zero. Initially we have two stable critical points, one of 
which is a minima and the other a maxima, and two hyperbolic points. As we increase h 2 a pair of critical points, one 
stable and the other hyperbolic, appears in the region around the maxima. The separatrix passing through this new 
hyperbolic point encloses both the maxima and the new stable point which is also a maxima. The two maxima may 
not be degenerate but clearly, for a range of energy values greater than the separatrix energy and less than the energy 
of the smaller maxima, we will have disconnected energy contours. Since the extra pair of critical points emerge from 
a single point, the phase diagram can be obtained by finding the condition for two roots of eqn (6) to be identical. 
The equation that we thus get can be solved numerically and we obtain the phase diagram shown in fig. 5. This 
shows various constant <f> sections of the phase diagram. In fig. 6 we plot constant energy contours for non zero <fi 
and for field values on the two sides of the phase boundary. The plot illustrates the appearance of critical points and 
associated with it the disconnected energy contours. 

To summarize the results obtained so far, we have shown that for a two spin system, for generic values of external 
fields, the dynamics is crgodic on the geometrically connected parts of energy contours. There arc as many sectors as 
disconnected parts. Studying the critical points enables us to obtain full information on the number of sectors as a 
function of the fields and energy. 

IV. ERGODICITY IN AN N SPIN SYSTEM 

We now give some results that are valid for any number of spins. In this case too, if there are closed disconnected 
surfaces with the same energy, the dynamics does not allow transitions between them. This simply follows from the 
fact that any line, along which only 9k varies, will intersect a closed energy surface at two points. Since we allow 
single flips only, at every time step, one can move only between these two points on the same surface. 

It can also be shown that for the many spins case, phase space trajectories obtained for generic initial conditions 
densely fill the connected parts of the constant energy surface. This is an important result since it implies ergodicity 
within sectors. Consider any point P(6\, 82, ---On) on the energy surface. Suppose we take any two nearest neighbor 
spins 9k and 0i and flip them successively while keeping all other spins fixed. The fixed spins in this case act as 
external fields. We know from the two spin case that for almost every value of the external fields, which in this 
case means for almost any configuration of the fixed spins, the trajectory will be dense in the restricted phase space. 
The restricted phase space is the intersection between the energy surface and the 9k-9i plane. By choosing different 
nearest neighbor pairs of spins, we can have (N — 1) independent directions emerging from the point P on the (N — 1) 
dimensional energy surface, along each of which there exists a dense orbit. This in turn means that, for every point 
there exists a neighborhood such that all points of the neighborhood can be reached by the dynamics. Hence the 
trajectory obtained by evolving any initial condition cannot have a boundary and so has to fill all of the connected 
part of the energy surface. 

Finally to obtain information on the connectivity of the energy surfaces, one can study the critical points. However 
unlike the two spin case, for more number of spins, we are not able to do this in full detail. 

In the following subsections, first as an illustrative example, we treat the case of three spins in external fields. Then 
we consider spins on general hypercubic lattices. 

A. Ergodicity in Three-spin case 

For an open chain of three spins interacting with each other and with external fields, the energy function is given 
by 

H = -hicos(9i - 4>i) - h2Cos(9 2 ) — h 3 cos(8 3 — </> 3 ) - cos(9i - 9 2 ) - cos(9 2 - 9 3 ) 

In this case depending on the the sequence in which the spins are flipped we get different maps in the three 
dimensional phase space. Here we look at the stochastic map in which any one spin is chosen with equal probability 
and flipped. Thus if we represent the three spins by complex numbers z\, z 2 and z 3 then at every time step we choose 
one of the following maps: 
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^ ~ /He"* + 1/Z2 

, _ h 2 + zi + z 3 
" 2 - / i2 + l/z 1 + l/z 3 1/Z2 

-l/*3 (7) 



J ft 3 e-^3 + l/z 2 " 

Let us study this dynamics for the special case of fields given by: hi = h 2 = /13 = 1.0 and <pi = (f> 3 = 0.0. We first 
look at the nature of the constant energy surfaces. For this we study the critical points. 

The critical points are obtained by setting to zero the first derivatives of 7i w.r.t 61,62 and 83. This is equivalent 
to the condition that at every site the net field is either parallel to the spin at that site or it vanishes. Thus from (7) 
we get: 

2 1 + z 2 

Zl = 



7 2 _ 
z 2 — 



1 + 1/Z2 

1 + zi + z 3 

f 1/zi + 1/2:3 

2 = 1 + Z 2 

3 1 + 1/Z2 ' 



(8) 



Solving the above equations we obtain the following critical points: Pi(l, 1, 1), P 2 (— 1, 1, —1), p3(L 1, —1), Pi(— 1, 1, 1), 
Ph{z\,—1, l/zi), Pe(zi, — 1, — zi). In fig. 7, we show the locations of these critical points in the (61,82, 63) phase space. 

By expanding Ji around the critical points it can be shown that Pi is a minimum of the energy with e = — 5 while 
P2 is a maximum with e = 3. The points P3 and P4 are hyperbolic points with e = — 1. The lines P5 and Pq consist 
of hyperbolic points and correspond to e = 1. Knowing the nature of all the critical points and by carefully looking at 
the energy function, we construct the following picture of the energy surfaces: The point Pi is the lowest energy state 
corresponding to all spins aligned along the field direction. For e w —5, the constant energy surfaces are ellipsoidal, 
enclosing the point Pi. In fact for —5 < e < —1 the energy surfaces are all closed deformed ellipses. At e = —1 the 
energy surface passes through the two unstable points P3 and P4. For — 1 < e < 1 we get spanning surfaces. The 
e = l surface is again closed and passes through the critical lines P5 and Pg. Finally for 1 < e < 3, the surfaces are 
closed around the other stable point at P 2 . We note that in all cases the energy surfaces are connected. 

Let us now look at the numerically obtained phase space trajectories for different energy values. These are shown 
in fig. 8. For every energy value the trajectory was obtained by evolving a single initial condition. In all cases the 
dynamics appears to be ergodic. The energy values were chosen so that all the different types of surfaces discussed 
above are generated. For example it can be seen that at e = —0.832 we get a spanning surface while for e > 1 the 
surfaces are closed. The surface in fig. 8c is a spanning surface topologically equivalent to the one in fig. 8d. To 
further verify the ergodic nature of the trajectories we plot in fig. 9 constant 82 sections of the phase space trajectory 
for the energy value e = —0.832. It is easy to understand the form of the contours if we rewrite the energy in the 
following way 

H = 2 cos(6» 2 /2) (003(6!! - 2 /2) + cos(6» 3 - 6> 2 /2)) + cos(6> 2 ). 

For other values of energy we get results similar to those in fig. 9. 

Thus in this three spin system we find that all energy surfaces are connected and the phase space trajectories 
densely fill the energy surfaces. Hence we have complete ergodicity. 

B. Spins on a lattice 

On a general hypercubic lattice consisting of N sites the equations for the critical points is 



- >Zj l/zi i = l,2 N. (9) 



1/ Zj 



The sum over j is over nearest neighbors of i. In all cases we assume periodic boundary conditions. We have not 
been able to obtain the full set of solutions of these equations and analyze their structure. We have done this for a 
class of solutions and this is described below. The one dimensional case is a special case and we first describe it. 
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(i) One dimensional lattice 

First we note that in this case, if we set all the coupling constants strictly to one then the dynamics is not ergodic 
even within sectors. The proof of ergodicity in III breaks down. This is because fixing all spins except two nearest 
neighbors now corresponds to the special case h\ = h 2 = 1.0 mentioned in section III, which was shown to be 
periodic. However it was seen that if we make the couplings different from one by arbitrarily small amounts, then we 
get ergodicity (within sectors). Assuming this has been done we can still determine the phase space structure with 
the couplings set equal to one. 

The set of eqns. (9) can be written in the following equivalent form: 

zl = z k - 1 z k+1 or z k _i + z k+ i = k=l,2...N. (10) 

Due to periodic boundary condition we have zq — zn and z^v+i = z\. Thus we have 2 N sets of equations whose 
solution would give us the critical points. It is easy to see that the only stable solutions (stable in all directions except 
for a uniform rotation) are obtained from the set of equations 

z 2 k = z k - 1 z k+1 fc = l,2.JV. (11) 

Every other of the 2^ sets gives saddle points since there are one or more sites where the field vanishes and at those 
sites one can freely rotate the spins without changing the energy. For this set since the fields don't vanish at any site 
they also correspond to fixed points of the dynamics. That the above set is stable is proved by explicitly solving for 
the fixed points and analyzing their nature. We get N solutions given by z k (n) — e 27Tlkn / N Zl with n going from 1 
to N. The factor z\ signifies the fact that the direction of one of the spins can be fixed arbitrarily. We now look at 
the matrix M k i = d 2 Ti./dd k ddi evaluated at the fixed points. The non-vanishing matrix elements thus obtained are 
M kk (n) = —2cos(2nn/N), M kk _i = M kk+ i = cos(2nn/N). Here n refers to the nth fixed point. This can be easily 
diagonalised and gives the eigenvalues 

Xg(n) = -2cos(27rn/A0(l - cos(2nq/N)) q = 0,l, —N - 1. 

Except the first eigenvalue q = which corresponds to a uniform rotation, all other eigenvalues are non- vanishing and 
have the same signs. Thus all fixed points are stable and the energy surfaces around them are elliptical. For the case 
when TV is a multiple of 4, a special case arises for the n = (7V/4)th fixed point. In this case all eigenvalues vanish. 
This is because in this case the field at every site vanishes and spins can rotate freely. 

The energy of the nth fixed point is given by e = TV cos(2im/N) and we note that every fixed point except those 
corresponding to the lowest and highest energies are doubly degenerate (for even TV). The energy surfaces around 
these fixed points thus consist of two disconnected parts each of which is closed. 

Another set of solutions of eqns. (10) is obtained by choosing any arbitrary direction and aligning every spin either 
along or opposite to it. The all up and the alternating up and down configurations have already been generated in 
the previous set of solutions. All the other solutions are saddles and it is easy to see that the degenerate points are 
all connected. The degeneracy of each solution is given by the number of configurations with a fixed number of bad 
bonds. 

In summary, for the one dimensional lattice, we have obtained all the fixed points and shown that they are stable 
and doubly degenerate. 



(ii) General case of d-dimensional hypercubic lattice of length L 

In higher dimension, eqns. (9) are highly nonlinear and classifying all solutions is a difficult task. However we 
obtain by inspection a set of solutions corresponding to spin-waves. These are of the form 

z k (n) = e 2nfk ^ L z (12) 

where the vector k denotes points on the d-dimensional lattice with integral coordinates (hi, ■■■kd) and n = 
(m, n,2, ■■■rid) denotes one of the L d solutions. The indices fcj and rn take values 1,2,...L. As before the solutions 
contain an arbitrary constant phase factor z. 

To understand the nature of the critical points we again look at the eigenvalues of the matrix M k i — d 2 TL/ dO^dOj, 
evaluated at the critical points. The non-vanishing matrix elements are in general given by M kk = — J2le<k> cos (% ~~ 
6j) and M kkp — cos(9 k — &k p )- Here k p denotes the nearest neighbors along the pth axis, namely the points 

(fei, &2, ■■■■kp ± 1, ...kd). In our case the matrix elements at the critical point n are M kk (h) = —2 J2 P =i co$,(2irn p / L) 
and M kk (n) = cos(2irn p / L) . It is easy to show that the eigenvalues of this matrix are 
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d 

\g(n) = -V]2cos(27rnp/L)(l - cos(2irq p /L)) q p = 0, ...N — 1. (13) 
P =i 

Unlike in the one-dimensional case eigenvalues corresponding to a critical point can have opposite signs. In fact the 
only stable critical points are those for which cos(2Tm p / L) has the same sign for all p. To show this suppose there 
exists in the set {n\,ri2---nd) two values n a and rib for which cos(2nn a /L) and cos(27rnfc/L) have the opposite signs. 
Then the eigenvalues for q = (0, 0, ...q a , 0, ..0) and q — (0, 0, qb, 0, ...0) have opposite signs. To find whether the 
critical points correspond to fixed points we look at the fields at each site. This is given by H^(n) = z k p (n) — 

2e 2m ^ n ^2 p : =1 cos(2irn p /L). Thus except for the small subset of points satisfying ^2 p=1 cos{2im p /L) = 0, all other 
critical points are fixed points. Unlike the one dimensional case, here we may have other fixed points corresponding 
to other solutions of eqns.(9). 

The energy of the critical point n is e = JV^ =1 cos(2irn p /L). We note that the spectrum of energies spans the 
entire allowed energy range. If d g is the number of non-zero coordinate values in the set n, then there is a degeneracy 
of 2 d <> associated with the corresponding energy value since we can flip the signs of rij without affecting the energies. 

The fact that there exist degenerate stable fixed points implies that there are disconnected energy surfaces corre- 
sponding to the neighborhoods of these points. This proves that the dynamics is non-ergodic for certain energies. It 
is however possible that we can get ergodicity for energies corresponding to unstable points and their neighborhood. 
The following numerical simulations address this question. 

We took a 20 x 20 square lattice and computed 1st, 2nd, 3rd and 4th nearest neighbor correlation functions 
(c(l),c(2),c(3),c(4)) using monte-carlo simulation. We study the two cases when the initial data is chosen close to a 
stable or close to an unstable point. In each case we compare the results with those obtained from randomly chosen 
initial configuration. Close to a stable fixed point the correlation function < cos(% — 0f. + fi) > can be calculated 
exactly. At the fixed point n this is simply cos(2ir R.n / L). The correlation c(r) evaluated in the simulations corresponds 
to (cos(2irrni) + cos(27rrn 2 ))/2. The following table gives the numerical data and also the exact values at the stable 
points: 



Initial Condition 


6/2 = c(l) 


c(2) 


c(3) 


c(4) 


n « (2,6)/Unstable 


0.24855 


0.06646 


0.01861 


0.00524 


Random 


0.24855 


0.06647 


0.01898 


0.00574 


n « (3,4)/Stable 


0.44640 


-0.55508 


-0.87305 


-0.24811 


Random 


0.44640 


0.23264 


0.13246 


0.07996 


Exact 


0.44840 


-0.55902 


-0.88004 


-0.25 



As expected for the stable case we get results totally different from the random case. On the other hand at the 
unstable points we get the same correlations as with random initial conditions. Thus as far as correlation functions 
are concerned there seems to be no dependence on initial conditions provided we don't start from near stable points. 

Next we show that the correlation functions obtained with randomly chosen initial conditions are the same as those 
obtained from a canonical simulation. In fig. 10 we plot the correlations c(2), c(3) and c(4) against the average 
energy given by c(l) for data obtained from the energy conserving dynamics and from a canonical simulation. In the 
canonical simulation each spin is coupled to a heat bath at some temperature and clearly we have ergodicity. For 
the energy conserving dynamics we started with a high energy configuration of spins pointing in random directions. 
Subsequently the energy was decreased by choosing spins randomly and aligning them along the local field direction. 
For each configuration obtained in this way we calculated the correlation functions. The data from the two approaches 
seem to match closely. 

The above numerical results suggest that: 

(1) the set of stable points and the sectors associated with them seem to have a small measure. 

(2) the unstable points with the same energies are either all connected by the dynamics or correlation functions in 
each sector are the same. 

A more complete study of the phase space structure is needed in order to understand the above results. 

Finally we note that as in the one dimensional case, another class of solutions of the critical point equations can be 
written. These are the solutions for which every spin points either along or opposite to a fixed direction. These critical 
points may be either stable or unstable. For example consider the spin configuration of a single up-spin in a sea of 
down-spins. This is an unstable fixed point. The different fixed points, corresponding to the N different locations of 
the single up-spin, are all connected. To show this suppose we start from a configuration which is arbitrarily close to 
such a fixed point. Now we fix all spins except for the single up-spin and any one of its down neighbors. The results 
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in section III imply that by flipping these two spins successively we can reach the configuration in which the two spins 
have interchanged their orientations. Repeating this for different nearest neighbor spin pairs, we can reach all the N 
fixed points. This is different from the Ising model with similar dynamics where the system would remain stuck in 
one of the fixed points. The above procedure can be used to prove hyperbolicity of other fixed points also. 

In conclusion, we have studied a stochastic energy conserving single spin flip dynamics of the XY model. We find 
that phase space trajectories fill connected parts of the energy surface and there is loss of crgodicity whenever the 
energy surface consists of closed disconnected parts. We also show that one can obtain information about the nature 
of the energy surfaces by studying the critical points of the energy function. This can be done in full detail for a 
two-spin system. For spins on hypercubic lattices we have proved that the dynamics is not ergodic for all energies and 
time averages of physical quantities depend on choice of initial conditions. However, the total weight of non-typical 
initial conditions appears to be small and the dynamics is effectively ergodic. 

I would like to thank Deepak Dhar for valuable suggestions all through the course of this work. 
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Figure Captions 

FIG. 1. Schematic picture of constant energy contours for two spin case. 



FIG. 2. Phase space trajectories of the map (3) for hi — 1, hi = 1.001, and 4> — 0. 

FIG. 3. Phase diagram for (j> — 0, showing nature of critical points for different field values. In regions Ri, R2 and -R3 there 
are four critical points of which two are hyperbolic, one a maxima and the other a minima. In region R4 there are six critical 
points. Three of them are hyperbolic, two are maxima and one is a minima. 

FIG. 4. Constant energy contours obtained by iteration of the map (3) starting from different initial conditions. The figures 
(a)-(e) are for five different field values corresponding to the points A-E in fig. 3. In all cases, some of the initial conditions 
were chosen close to the hyperbolic points, so that the separatrices can be seen. 

FIG. 5. Phase diagram for different values of <j>. In every case, for field values in the region bounded by the three curves, 
there is a range of energy where there are two sectors. 

FIG. 6. Constant energy contours for <j> = 7r/16. The field strengths hi and hi were chosen to lie on two sides of the phase 
boundary shown in fig. 5. In (a) hi = 1.25, hi = 0.77 while in (b) hi — 1.25, /12 = 0.89. 

FIG. 7. Position of critical points of three spin energy function in the (#i,#2,#3) phase space. The value of 9 ranges from 
— 7r to 7r along each of the axes. 
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FIG. 8. Three dimensional phase space trajectories for field values hi = h 2 = /13 = 1.0. Each of the 5 figures was obtained 
by iterating different initial conditions corresponding to different energy values. 



FIG. 9. Constant 2 sections of the phase space trajectories for the energy value e = —0.832. 



FIG. 10. Correlation functions c(2),c(3) and c(4) plotted against c(l) for data obtained using energy conserving and 
non-conserving dynamics. In each case the averages were taken over 10 5 monte-carlo steps. 
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